library("msy")

# Blim = 0.2 * SSBmax
result.brp.blim <- max(ssb(result.xsa.stock)) * 0.2
# get fbar values as vector
fvals <- as.vector(fbar(result.xsa.stock))

# found extremes of whiskers by boxplot status to remove unreal data
wstats <- boxplot.stats(fvals)$stats
# remove values out of box of extremums 
fvals <- fvals[fvals > wstats[[1]]]
fvals <- fvals[fvals < wstats[[5]]]

# get variance (standard deviation) of Fbar
fsd <- sd(fvals)

# Bpa = Blim * e^(1.645*stdev)
result.brp.bpa <- result.brp.blim * exp(qnorm(0.95)*fsd)


# fit SR models via eqsim
result.eqsim.srmodel <- eqsr_fit(result.xsa.stock, 
                                 n = 5000, 
                                 models = c("Ricker", "Bevholt", "Segreg")
                                 #remove.years = c(2013,2014,2015)
                                 )

lastYear <- as.numeric(result.xsa.stock@range["maxyear"])

# calculate reference point based on model
result.eqsim.brp <- eqsim_run(result.eqsim.srmodel, 
                              bio.years = c(lastYear-3, lastYear), 
                              extreme.trim = c(0.05, 0.95), 
                              Bpa = result.brp.bpa,
                              Blim = result.brp.blim)

# create new FLR sr-model, segreg from eqsim params
tmp.model <- filter(result.eqsim.srmodel$sr.det, model=="Segreg")
result.sr.segreg <- FLSR(model = "segreg")
result.sr.segreg@params["a"] <- as.numeric(tmp.model["a"])
result.sr.segreg@params["b"] <- as.numeric(tmp.model["b"])

# get mean recruitment of latest 5 years without terminal year
tmp.rec.mean <- exp(mean(log(window(rec(result.xsa.stock), 2014, 2018))))
# build geomean SR model
result.sr.geomean <- FLSR(model = "geomean", params=FLPar(tmp.rec.mean))

# build test model
test.geomean <- fwd(data.stf, sr=result.sr.geomean, target.forecast.F_SQ)
test.segreg <- fwd(data.stf, sr=result.sr.segreg, target.forecast.F_SQ)
